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Abstract 

We  present  a  new  algorithm  to  model  the  input  uncertainty  and  its  propagation  in  incompressible  flow  simu¬ 
lations.  The  stochastic  input  is  represented  spectrally  by  employing  orthogonal  polynomial  functionals  from  the 
Askey  scheme  as  trial  basis  to  represent  the  random  space.  A  standard  Galerkin  projection  is  applied  in  the 
random  dimension  to  obtain  the  equations  in  the  weak  form.  The  resulting  system  of  deterministic  equations  is 
then  solved  with  standard  methods  to  obtain  the  solution  for  each  random  mode.  This  approach  can  be  consid¬ 
ered  as  a  generalization  of  the  original  polynomial  chaos  expansion,  first  introduced  by  N.  Wiener  (1938).  The 
original  method  employs  the  Hermite  polynomials  (one  of  the  thirteen  members  of  the  Askey  scheme)  as  the 
basis  in  random  space.  The  algorithm  is  applied  to  pressure-driven  channel  flows  with  random  wall  boundary 
conditions,  and  to  external  flows  with  random  freestream.  Efficiency  and  convergence  are  studied  by  comparing 
with  exact  solutions  as  well  as  numerical  solutions  obtained  by  Monte  Carlo  simulations.  It  is  shown  that  the 
generalized  polynomial  chaos  method  promises  a  substantial  speed-up  compared  with  the  Monte  Carlo  method. 
The  utilization  of  different  type  orthogonal  polynomials  from  the  Askey  scheme  also  provides  a  more  efficient  way 
to  represent  general  non-Gaussian  processes  compared  with  the  original  Wiener-Hermite  expansions. 
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1  Introduction 


Recently  there  has  been  an  intense  interest  in  verification  and  validation  of  large-scale  simulations  and  in  modeling 
uncertainty  [1,  2,  3].  In  simulations,  just  like  in  the  experiments,  we  often  question  the  accuracy  of  the  results  and 
we  construct  a  posteriori  error  bounds,  but  the  new  objective  is  to  model  the  uncertainty  from  the  beginning  of 
the  simulations  and  not  simply  as  an  afterthought.  Numerical  accuracy  and  error  control  have  been  employed  in 
simulations  for  some  time  now,  at  least  for  the  modern  discretizations,  e.g.  [4,  5].  However,  there  are  always  some 
uncertain  components  associated  with  the  physical  problems,  specifically  with  such  diverse  factors  as  constitutive 
laws,  boundary  and  initial  conditions,  transport  coefficients,  source  and  interaction  terms,  geometric  irregularities 
(e.g.  roughness),  etc.  . 

Most  of  the  research  efforts  in  CFD  research  so  far  have  been  in  developing  efficient  algorithms  for  different 
applications,  assuming  ideal  inputs  with  precisely  defined  computational  domains.  With  the  field  reaching  some 
degree  of  maturity  now,  we  naturally  pose  the  more  general  question  of  how  to  model  uncertainty  and  stochastic 
inputs,  and  how  to  formulate  algorithms  to  accurately  reflect  the  propagation  of  the  uncertainty.  To  this  end,  the 
Monte  Carlo  approach  can  be  employed  but  it  is  computationally  expensive  and  is  only  used  as  the  last  resort. 
The  sensitivity  method  is  a  more  economical  approach,  based  on  the  moments  of  samples,  but  it  is  less  robust  and 
depends  strongly  on  the  modeling  assumptions  [6].  One  popular  technique  is  the  perturbation  method  where  all 
the  stochastic  quantities  are  expanded  around  their  mean  via  Taylor  series.  This  approach,  however,  is  limited  to 
small  perturbations  and  does  not  readily  provide  information  on  high-order  statistics  of  the  response.  The  resulting 
system  of  equations  becomes  extremely  complicated  beyond  second-order  expansion.  Another  approach  is  based  on 
expanding  the  inverse  of  the  stochastic  operator  in  a  Neumann  series,  but  this  too  is  limited  to  small  fluctuations, 
and  even  combinations  with  the  Monte  Carlo  method  seem  to  result  in  computationally  prohibitive  algorithms  for 
complex  systems  [7]. 

A  more  effective  approach  pioneered  by  Ghanem  and  Spanos  [8]  in  the  context  of  finite  elements  for  solid  mechanics 
is  based  on  a  spectral  representation  of  the  uncertainty.  This  allows  high-order  representation,  not  just  first-order  as 
in  most  perturbation-based  methods,  at  high  computational  efficiency.  It  is  based  on  the  original  theory  of  Wiener 
(1938)  on  homogeneous  chaos  [9,  10].  This  approach  was  employed  in  turbulence  in  the  1960s  [11,  12,  13].  However, 
it  was  realized  that  the  chaos  expansion  converges  slowly  for  turbulent  field  [14,  15,  16],  so  polynomial  chaos  did  not 
receive  much  attention  for  a  long  time. 

The  main  purpose  of  this  paper  is  to  demonstrate  that  the  polynomial  chaos  expansion  can  be  effective  in  modeling 
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uncertainties  associated  with  fluid  flows.  When  the  uncertainty  has  relatively  strong  correlation,  the  chaos  expansion 
converges  fast;  in  the  ideal  case  it  converges  exponentially  fast  due  to  the  fact  that  it  is  a  spectral  expansion  in  the 
random  space.  The  spectral  representation  of  the  uncertainty  is  based  on  a  trial  basis  where  9  denotes 

the  random  event.  For  example,  the  vorticity  has  the  following  finite-dimensional  representation 

p 

w(x,i,0)  =  ^^(x,*)^^#)). 

»= 0 

Here  w,  (x,  t)  represents  the  deteministic  coefficents  and  will  be  denoted  as  the  random  mode  (?)  of  the  vorticity.  The 
random  trial  basis  is  a  set  of  complete  orthogonal  polynomials  in  terms  of  the  multi-dimensional  random  variable 
£(0)  with  a  specific  probability  distribution.  \fr(£(0))  is  a  functional,  as  it  is  a  function  of  random  variables  £  which 
are  functions  of  the  random  parameter  9  €  [0,1].  For  the  original  polynomial  chaos  introduced  by  Wiener,  the 
polynomial  trial  basis  is  the  Hermite  polynomials  in  terms  of  multi-dimensional  Gaussian  random  variables.  In  this 
paper,  we  will  apply  this  expansion  to  fluid  flows  and  further  generalize  the  trial  basis  to  other  orthogonal  polynomials 
from  the  Askey  scheme  [17].  For  different  types  of  basis  polynomials,  the  random  variables  £(0)  are  not  restricted 
to  the  Gaussian  variables.  Therefore,  we  have  additional  flexibility  to  represent  the  non-Gaussian  processes  more 
efficiently.  The  theory  of  orthogonal  functionals  plays  a  key  role  in  the  algorithms  developed  here.  We  note  that  the 
Monte  Carlo  algorithm  can  be  thought  of  as  a  subcase  of  the  above  representation  corresponding  to  the  collocation 
procedure  where  the  test  basis  is  \Eq(0)  =  6(9  —  9i ),  where  6  is  the  Kronecker  delta  function  and  9t  refers  to  an 
isolated  random  event. 

The  algorithms  we  develop  here  are  general  but  we  present  applications  with  uncertainty  associated  with  boundary 
conditions.  This  situation  is  encountered,  for  example,  in  micro-channel  flows  but  also  in  classical  flows  such  as  the 
freestream  flow  past  bluff  bodies.  The  generalized  polynomial  chaos  expansion  can  handle  both  Gaussian  and 
non-Gaussian  random  processes.  For  certain  distributions  there  exist  a  “best”  representation  which  results  in  fast 
convergence  rate.  For  example,  for  Poisson  distributions  is  the  Charlier  polynomials,  for  Gamma  distributions  the 
Laugerre  polynomials,  for  binomial  distributions  the  Krawtchouk  polynomials,  for  the  beta  distributions  the  Jacobi 
polynomials,  etc.  . 

In  the  next  section  we  review  the  theory  of  the  Askey  scheme  of  hypergeometric  orthogonal  polynomials,  and  in 
section  3  we  present  the  framework  of  the  generalized  polynomial  chaos.  In  section  4  we  address  its  implementation 
details  when  applied  to  Navier-Stokes  equations.  In  section  5  we  present  the  computational  results  of  various 
applications  and  demonstrate  the  convergence  property  of  the  chaos  expansion.  We  conclude  the  paper  with  a 
discussion  on  open  questions.  In  the  appendix  we  include  a  brief  review  of  orthogonal  polynomials. 
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2  The  Askey  Scheme  of  Hypergeometric  Orthogonal  Polynomials 


The  theory  of  orthogonal  polynomials  is  relatively  mature  and  several  books  have  been  devoted  to  their  study  (e.g., 
[18,  19,  20]).  More  recent  work  has  shown  that  an  important  class  of  orthogonal  polynomials  belongs  to  the  Askey 
scheme  of  the  hypergeometric  polynomials  [17].  In  this  section,  we  briefly  review  the  theory  of  hypergeometric 
orthogonal  polynomials;  we  adopt  the  notation  of  [21,  22]. 

2.1  The  Generalized  Hypergeometric  Series 


We  first  introduce  the  Pochhammer  symbol  (a)n  defined  by 

1,  if  n  =  0, 


(@)n  — 


In  terms  of  Gamma  function,  we  have 


r (a)  ’ 

The  generalized  hypergeometric  series  rFs  is  defined  by 


i(a  +  1)  •  •  •  (a  +  n  —  1),  if  n  =  1, 2,  3, . . . . 


T(a  +  n) 

( a)n  =  , — ,  n  >  0. 


jp  f  ,  ,  n  _  (ai)fc  ’  ’ '  (ar)k  z 

’  r’  l!"'  ’bs'z)~2^  (h)k---(bs)k  k\ 


k= 0 


(i) 


(2) 


(3) 


where  bi  y  0,  —1,  —2, . . .  for  i  =  {1, . . . ,  s}  to  ensure  the  denominator  factors  in  the  terms  of  the  series  are  never 
zero.  The  radius  of  convergence  p  of  the  hypergeometric  series  is 


oo  if  r  <  s  +  1, 
p  =  {  1  if  r  =  s  +  1, 

0  if  r  >  s  +  1. 


(4) 


Some  elementary  cases  of  the  hypergeometric  series  are:  the  exponential  series  qFq  and  the  binomial  series  iFq. 

If  one  of  the  numerator  parameters  a*,  i  =  1, . . . ,  r  is  a  negative  integer,  say  ci\  =  —n,  the  hypergeometric  series 
(3)  terminates  at  the  nth- term  and  becomes  a  polynomial  in  z, 

( -n)k  •  •  •  ( ar)k  zk 


rFs (  n,  *  *  •  ,  ctr ,  ,  ,bs,  z)  —  ^  ] 


fc= o 


(bi)k- ■  ■  {bs)k  k\' 


(5) 


2.2  Some  Properties  of  the  Orthogonal  Polynomials 


A  system  of  polynomials  {Qn{x),  n  £  J\f}  where  Qn(x )  is  a  polynomial  of  exact  degree  n  and  J\f  =  {0,1,2,...}  or 
A f  =  {0, 1, . . . ,  IV}  for  a  finite  nonnegative  integer  N,  is  an  orthogonal  system  of  polynomials  with  respect  to  some 
real  positive  measure  </>  if  the  following  orthogonality  relations  are  satisfied 

I  Qn(x')Qm{^x)d(j)(^x')  —  hn6nrn ,  77.,  777.  £  AJ" ,  (6) 

Js 
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where  S  is  the  support  of  the  measure  <j>  and  the  hn  are  non-zero  constants.  The  system  is  called  orthonormal  if 
hn  —  1. 

The  measure  <f>  often  has  a  density  w(x)  or  weights  w(i)  at  points  aq  in  the  discrete  case.  The  relations  (6)  then 
become 

/  Qn{x)Qm{x)w(x)dx  =  h2n5nrn,  71,  TO  S  Af ,  (7) 

Js 

in  the  continuous  case,  or 

M 

Qn(xi)Qm(xi)w(xi)  =  h2n8nm,  n,  m  G  A/*,  (8) 

i= 0 

in  the  discrete  case  where  it  is  possible  that  M  =  oo. 

The  density  w(x),  or  weights  w(i)  in  the  discrete  case,  is  also  commonly  referred  as  the  weighting  function  in  the 
theory  of  orthogonal  polynomials.  It  will  be  shown  later  that  the  weighting  functions  for  some  orthogonal  polynomials 
are  identical  to  certain  probability  functions.  For  example,  the  weighting  function  for  the  Hermite  polynomials  is 
the  same  as  probability  density  function  of  the  Gaussian  random  variables.  This  fact  plays  an  important  role  in 
representing  stochastic  processes  with  orthogonal  polynomials. 

2.3  The  Askey  Scheme 


The  Askey  scheme,  which  is  represented  as  a  tree  structure  in  figure  1  (following  [22]),  classifies  the  hypergeometric 
orthogonal  polynomials  and  indicates  the  limit  relations  between  them.  The  ‘tree’  starts  with  the  Wilson  polynomials 
and  the  Racah  polynomials  on  the  top.  They  both  belong  to  the  class  4F3  of  the  hypergeometric  orthogonal 
polynomials  given  by  equation  (5).  The  Wilson  polynomials  are  continuous  while  the  Racah  polynomials  are  discrete. 
The  lines  connecting  different  polynomials  denote  the  limit  transition  relationships  between  them;  this  implies  that 
the  polynomials  at  the  lower  end  of  the  lines  can  be  obtained  by  taking  the  limit  of  one  of  the  parameters  from  their 
counterparts  on  the  upper  end.  For  example,  the  limit  relation  between  Jacobi  polynomials  Pha’^\x)  and  Hermite 
polynomials  Hn(x)  is 


lim  a~^nPia’a) 


Hn(x) 
2  nn\ 


and  between  Meixner  polynomials  Mn(x:  ft,  c)  and  Charlier  polynomials  Cn(x',a)  is 


lim  Mn  (  x;  (3,  °  )  =  Cn(x;  a). 

p— ►  00  \  a  p  ) 

For  a  detailed  account  of  the  limit  relations  of  Askey  scheme,  the  interested  reader  should  consult  [21]  and  [22]. 

The  orthogonal  polynomials  associated  with  the  generalized  polynomial  chaos,  which  will  also  be  called  the 
Askey-Chaos  hereafter,  include:  Hermite,  Laguerre,  Jacobi,  Charlier,  Meixner,  Krawtchouk  and  Hahn  polynomials. 
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Figure  1:  The  Askey  scheme  of  orthogonal  polynomials 
A  review  of  their  definitions  and  properties  can  be  found  in  the  appendix  of  this  paper. 

3  The  Generalized  Polynomial  Chaos 

In  this  section  we  introduce  the  generalized  polynomial  chaos  expansion  along  with  the  Karhunen-Loeve  (KL)  ex¬ 
pansion,  another  classical  technique  for  representing  random  processes.  The  KL  expansion  can  be  used  in  some  cases 
to  represent  efficiently  the  known  stochastic  fields,  i.e.,  the  stochastic  inputs. 

3.1  The  Original  Wiener  Polynomial  Chaos:  Hermite-Chaos 

The  original  polynomial  chaos,  also  termed  as  the  homogeneous  chaos,  was  first  introduced  by  Wiener  [9].  It  employs 
the  Hermite  polynomials  in  terms  of  Gaussian  random  variables.  According  to  a  theorem  by  Cameron  &  Martin 
[23],  it  can  approximate  any  functionals  in  ^(C)  and  converges  in  the  L2(C)  sense,  where  C  is  the  space  of  real 
functions  which  are  continuous  on  the  interval  [0, 1]  and  vanish  at  0.  Therefore,  polynomial  chaos  provides  a  means 
for  expanding  second-order  random  processes  in  terms  of  Hermite  polynomials.  Second-order  random  processes  are 
processes  with  finite  variance,  and  this  applies  to  most  physical  processes.  Thus,  a  general  second-order  random 
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process  X(9),  viewed  as  a  function  of  6 ,  i.e.  the  random  event,  can  be  represented  in  the  form 


X{8)  =  a0H0 

OO 

+  Y  Qii#l(&i(0)) 

*i=l 
oo  i\ 

+  EE  -^2  (&  (*),&(*)) 

il=l  22  =  1 

OO  2l  22 

+  EE  E«u^3^3teiW,^w,^3W) 

*1=1  *2  =  1  *3  =  1 

+  •••,  (9) 


where  , . . . ,  ^,n)  denote  are  Hermite  polynomials  of  order  n  in  terms  of  the  multi-dimensional  independent 

standard  Gaussian  random  variables  i  =  (£jx , . . . ,  £,n)  with  zero  mean  and  unit  variance.  The  above  equation  is  the 
discrete  version  of  the  original  Wiener  polynomial  chaos  expansion,  where  the  continuous  integrals  are  replaced  by 
summations.  The  general  expression  of  the  Hermite  polynomials  is  given  by 


5  *  *  •  5  £*n) 


(-i  r 


dn 


5£u  '  '  '  diin 


For  example,  the  one-dimensional  Hermite  polynomials  are: 


(10) 


*0  =  1,  *1=£,  *2  =  £2-l,  *3  =  £3-3£, 


(11) 


For  notational  convenience,  equation  (9)  can  be  rewritten  as 

OO 

*(*)  =  E>j*i(*),  (12) 

3=  0 

where  there  is  a  one-to-one  correspondence  between  the  functions  Hn(i ,...,£»„)  and  and  also  between  the 

coefficients  aj  and  Oj  jr.  In  equation  (9)  the  summation  is  carried  out  according  to  the  order  of  the  Hermite 
polynomials,  while  in  equation  (12)  it  is  simply  a  re-numbering  with  the  polynomials  of  lower  order  counted  first. 
For  clarity,  the  two-dimensional  expansion  is  shown  here,  both  in  the  fully  expanded  form  (  see  equation  (9)) 


X(9)  —  aoHo  +  aiHi(ii)  +  a2.f/i(£2) 

+  aii-ff2(£i)£i)  +  012^2(^2,  ^1)  +  «22  #2  (£2,  £2)  +  ■•• ,  (13) 


and  the  simplified  form  (see  equation  (12)) 


X(9)  =  aoTo  +  01*1  +  02*2  +  03*3  +  <23 ’F4  +  asTs  +  . . . 

=  a0  +  ai£i  +  a2^2  +  03 (£1  —  1)  +  a3(£i£2)  +  as(£2  ~  1)  +  ■  ■  ■  (14) 
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The  polynomial  chaos  forms  a  complete  orthogonal  basis  in  the  L2  space  of  the  Gaussian  random  variables,  i.e., 

<  TjTj  >=<  Tf  >  Sij,  (15) 

where  Sij  is  the  Kronecker  delta  and  <  •,  •  >  denotes  the  ensemble  average.  This  is  the  inner  product  in  the  Hilbert 
space  of  the  Gaussian  random  variables 

<  m9(t)  >=  j  mgiowiodz.  ae) 

The  weighting  function  is 

(17) 

where  n  is  the  dimension  of  What  distinguishes  the  Wiener-Hermite  expansion  from  many  other  possible  complete 
sets  of  expansions  is  that  the  polynomials  here  are  orthogonal  with  respect  to  the  weighting  function  W(£)  which 
has  the  form  of  the  multi-dimensional  independent  Gaussian  probability  distribution  with  unit  variance.  We  will  use 
the  term  Hermite-Chaos  hereafter  to  denote  the  Wiener  polynomial  chaos. 

3.2  The  Generalized  Polynomial  Chaos:  Askey-Chaos 

The  Hermite-Chaos  expansion  has  been  quite  effective  in  solving  stochastic  differential  equations  with  Gaussian 
inputs  as  well  as  certain  types  of  non-Gaussian  inputs,  e.g.  lognormal  distributions  [8],  [24,  25];  this  can  be  justified 
by  the  Cameron-Martin  theorem  [23].  However,  for  general  non-Gaussian  random  inputs,  the  convergence  rate  is  not 
fast.  In  some  cases  the  convergence  rate  is,  in  fact,  severely  deteriorated. 

In  order  to  deal  with  more  general  random  inputs,  we  introduce  the  generalized  polynomial  chaos  expansion, 
the  Askey-Chaos,  as  a  generalization  of  the  original  Wiener’s  Hermite-Chaos  expansion.  The  expansion  basis  of 
the  Askey-Chaos  is  formed  by  the  complete  set  of  orthogonal  polynomials  from  the  Askey  scheme  (see  section  2.3). 
Similar  to  section  3.1,  we  represent  the  general  second-order  random  process  X{6)  as 

X(9)  =  a^Io 

OO 

+  CuW) 

ii  =  l 
00  i\ 

ii=l  «2  — 1 
00  i\  12 

+  EEE  1^2  ^3 -^3  (C iMMO)M0)) 

*1=1  *2  =  1  *3=1 


(18) 


where  In(Ciir  ■  ■  •  >  0„)  denotes  the  Askey-Chaos  of  order  n  in  terms  of  the  multi-dimensional  random  variables  £  = 
(Cn)  •  •  ■  ?C*n)-  -^n  the  Askey-Chaos  expansion,  the  polynomials  In  are  not  restricted  to  Hermite  polynomials  but 
instead  they  could  be  any  member  of  the  Askey  scheme,  as  shown  in  figure  1.  Again  for  notational  convenience,  we 
rewrite  equation  (18)  as 

OO 

*(*)  =  £  WO.  (19) 

3=0 

where  there  is  a  one-to-one  correspondence  between  the  functions  In{Cin  ■  •  • ,  C*n)  and  'Lj  ( C) ,  and  their  coefficients 
Cj  and  Since  each  type  of  polynomials  from  the  Askey  scheme  form  a  complete  basis  in  the  Hilbert  space 

determined  by  their  corresponding  support,  we  can  expect  each  type  of  Askey-Chaos  to  converge  to  any  L2  functional 
in  the  L2  sense  in  the  corresponding  Hilbert  functional  space  as  a  generalized  result  of  Cameron-Martin  theorem 
([23]  and  [26]).  The  orthogonality  relation  of  the  Askey-Chaos  polynomial  chaos  takes  the  form 

<  >=<  >  Si:i ,  (20) 

where  <5^  is  the  Kronecker  delta  and  <  ■,  ■  >  denotes  the  ensemble  average  which  is  the  inner  product  in  the  Hilbert 
space  of  the  variables  ( 

<  f(C)g(C)  >=  J  mg«)W(C)dC,  (21) 

or 

<mg{Q>='52mg{QW{C)  (22) 

c 

in  the  discrete  case.  Here  W(C)  is  the  weighting  function  corresponding  to  the  Askey  polynomials  chaos  basis  {<!>;}; 
see  the  appendix  for  detailed  formulas.  Some  types  of  orthogonal  polynomials  from  the  Askey  scheme  have  weighting 
functions  of  the  same  form  as  the  probability  function  of  certain  types  of  random  distributions.  In  practice,  we  then 
choose  the  type  of  independent  variables  C,  in  the  polynomials  {$,(£)}  according  to  the  type  of  random  distributions 
as  shown  in  table  1.  It  is  clear  that  the  original  Wiener  polynomial  chaos  corresponds  to  the  Hermite-Chaos  and 
is  a  subset  of  the  Askey-Chaos.  The  Hermite-,  Laguerre-  and  Jacobi-Chaos  are  continuous  chaos,  while  Charlier-, 
Meixner-,  Krawtchouk-  and  Hahn-Chaos  are  discrete  chaos.  It  is  worth  mentioning  that  the  Legendre  polynomials, 
which  is  a  special  case  of  the  Jacobi  polynomials  (x)  with  parameters  a  =  (3  =  0,  correspond  to  an  important 

distribution  —  the  uniform  distribution.  Due  to  the  importance  of  the  uniform  distribution,  we  list  it  separately  in 
the  table  and  term  the  corresponding  chaos  expansion  as  the  Legendre-Chaos. 
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Random  inputs 

Wiener- Askey  chaos 

Support 

Continuous 

Gaussian 

Gamma 

Beta 

Uniform 

Hermite-Chaos 

Laguerre-Chaos 

Jacobi-Chaos 

Legendre-Chaos 

(  —  00,  00 ) 

[0,  oo) 

[a,b\ 

[a,  b\ 

Discrete 

Poisson 

Binomial 

Negative  Binomial 
Hypergeometric 

Charlier-Chaos 

Krawtchouk-Chaos 

Meixner-Chaos 

Hahn-Chaos 

{0,1,2,...} 
{0,1, •••A} 
{0,1,2,...} 
{0,1,- ..,iv} 

Table  1:  Correspondence  of  the  type  of  Wiener- Askey  polynomial  chaos  to  the  type  of  random  inputs  (N  >  0  is  a 
finite  integer). 


3.3  The  Karhunen-Loeve  Expansion 

The  Karhunen-Loeve  (KL)  expansion  [27]  is  another  way  of  representing  a  random  process.  It  is  based  on  the  spectral 
expansion  of  the  covariance  function  of  the  process.  Let  us  denote  the  process  by  /i(x,  9)  and  its  covariance  function 
by  Rhh(x,y),  where  x  and  y  are  the  spatial  or  temporal  coordinates.  By  definition,  the  covariance  function  is  real, 
symmetric,  and  positive  definite.  All  eigenfunctions  are  mutually  orthogonal  and  form  a  complete  set  spanning  the 
function  space  to  which  /i(x,  9)  belongs.  The  KL  expansion  then  takes  the  following  form: 

OO 

ft.(x,  9)  =  h(x)  +  ^2  V/A</>i(x)£l(6>),  (23) 

i—1 

where  /i,(x)  denotes  the  mean  of  the  random  process,  and  £i(9)  forms  a  set  of  independent  random  variables.  Also, 
0t(x)  and  Aj  are  the  eigenfunctions  and  eigenvalues  of  the  covariance  function,  respectively,  i.e. , 

jRhh(x,y)My)dy  =  XiMx).  (24) 

Among  many  possible  decompositions  of  a  random  process,  the  KL  expansion  is  optimal  in  the  sense  that  the 
mean-square  error  of  the  finite  representation  of  the  process  is  minimized.  Its  use,  however,  is  limited  as  the  covariance 
function  of  the  solution  process  is  often  not  known  a  priori.  Nevertheless,  the  KL  expansion  provides  an  effective 
means  of  representing  the  input  random  processes  when  the  covariance  structure  is  known. 

4  The  Askey-Chaos  for  Navier-Stokes  Equations 

In  this  section  we  present  the  solution  procedure  for  solving  the  stochastic  Navier-Stokes  equations  by  generalized 
polynomial  chaos  expansion.  The  randomness  in  the  solution  can  be  introduced  through  boundary  conditions,  initial 
conditions,  forcing,  etc.. 
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4.1  Governing  Equations 


We  employ  the  incompressible  Navier-Stokes  equations 


0,  (25) 

-Vn  +  t?e_1V2u,  (26) 

where  II  is  the  pressure  and  Re  the  Reynolds  number.  All  flow  quantities,  i.e.,  velocity  and  pressure  are  considered 
stochastic  processes.  A  random  dimension,  denoted  by  the  parameter  9 ,  is  introduced  in  addition  to  the  spatial- 
temporal  dimensions  (x,t),  thus 

u  =  u(x,  t;0);  II  =  n(x,f;0).  (27) 


V  •  u  = 


<9u 

9t 


+  (u  •  V)u  = 


We  then  apply  the  generalized  polynomial  chaos  expansion,  or  the  Askey-Chaos  (19),  to  these  quantities  and  obtain 


u(x,  t;  9)  =  ^Uj(x,£)$i(C(0));  n(x,t;0)  =  ^  11*  (x,  t)$i(C(0))> 


(28) 


i=0 


i— 0 


where  we  have  replaced  the  infinite  summation  in  infinite  dimension  of  C  hr  equation  (12)  by  a  truncated  finite-term 


summation  in  finite  dimensional  space  of  C-  The  total  number  of  expansion  terms,  ( P  +  1),  depends  on  the  number 
of  random  dimensions  (n)  of  C  and  the  highest  order  ( p )  of  the  polynomials  4*  [8]: 

p= yi  “i  ii(n+r)-  (29) 

s=l  '  r=0 

The  most  important  aspect  of  the  above  expansion  is  that  the  random  processes  have  been  decomposed  into  a  set 
of  deterministic  functions  in  the  spatial-temporal  variables  multiplied  by  the  random  basis  polynomials  which  are 
independent  of  these  variables. 

Substituting  (28)  into  Navier-Stokes  equations  ((25)  and  (26))  and  noting  that  the  partial  derivatives  are  taken 


in  physical  space  and  thus  commute  with  the  operations  in  random  space,  we  obtain  the  following  equations 

p 


u  =  o, 


i= 0 


^  9u4x,  *)  $,  +  ^  ^[(u. .  V)Ui)]$^  =  -^Vni(x,t)$i  +  i?e-1^V2ui$i. 


(30) 

(31) 


i= 0  ”  ”  z=0  j— 0  2=0  2=0 

We  then  project  the  above  equations  onto  the  random  space  spanned  by  the  basis  polynomials  {4>,}  by  taking  the 
inner  product  of  above  equation  with  each  basis.  By  taking  <  ■,  4>j.  >  and  utilizing  the  orthogonality  condition  (15), 


we  obtain  the  following  set  of  equations: 
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For  each  k  =  0, . . .  P, 


V-ufc  =  0,  (32) 

d  I  p  p 

eijfe[(uj  •  V)uj-)]  =  -VII k  +  Re  1V2ufe,  (33) 

k  i=0  j— 0 

where  e,^  =<  >.  Together  with  <  <Ff  >,  the  coefficients  can  be  evaluated  analytically  from  the 

definition  of  <&;.  The  set  of  equations  consists  of  (P  +  1)  system  of  ‘Navier-Stokes-like’  equations  for  each  random 
mode  coupled  through  the  convective  terms. 

4.2  Numerical  Formulation 

4.2.1  Temporal  Discretization 


We  employ  the  semi-implicit  high-order  fractional  step  method,  which  for  the  standard  deterministic  Navier-Stokes 
equations  ((25)  and  (26))  has  the  form  [28]: 


1  n-q 


u  —  V  n  a„u" 

9“°  -  =  -£/M(  u-v)u] 


At 


At 


9=o 


=  -vnn+1, 


7oun+1  —  u 

At 


=  Pe_1V2u"+1, 


(34) 

(35) 

(36) 


where  J  is  order  of  accuracy  in  time  and  a,  (3  and  7  are  integration  weights.  A  pressure  Poisson  equation  is  obtained 
by  enforcing  the  discrete  divergence-free  condition  V  •  u"+1  =  0 


V2nn+i  =  J_v  . 
At 


(37) 


with  the  appropriate  pressure  boundary  condition  given  as 


an 


—  =  -n  •  [u  +  Pe-1V  x  uP+1] 
on  L  J 


(38) 


where  n  is  the  outward  unit  normal  vector  and  u  =  V  x  u  is  the  vorticity.  The  method  is  stiffly-stable  and  achieves 
third-order  accuracy  in  time;  the  coefficients  for  the  integration  weights  can  be  found  in  [29]. 

In  order  to  discretize  the  stochastic  Navier-Stokes  equations,  we  apply  the  same  approach  to  the  coupled  set  of 
equations  (32)  and  (33): 
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For  each  k  =  0, . . . ,  P, 


Ufc  J2q=0  a9U, 

A t 


Eft 


<H>  n 
h  q— 0 


EE 


•i— 0  j= 0 


n—q 


Ufc  ~  tifc 
At 

7oufc+1  ~  u k 
At 


-vn”+1, 

1 Se_1V2<+1. 


(39) 

(40) 

(41) 


The  discrete  divergence-free  condition  for  each  mode  V  •  u^+1  =  0  results  in  a  set  of  consistent  Poisson  equations  for 
each  pressure  mode 

v2n"+1  =  Ev  •  ufc)  k  =  o,...,p,  (42) 

with  appropriate  pressure  boundary  condition  derived  similarly  as  in  [28] 


^  =  -n  •  [ufc  +  ite-'V  x  w”+1]  ,  k  =  0, . . .  ,P,  (43) 

where  n  is  the  outward  unit  normal  vector  along  the  boundary,  and  =  V  x  u k  is  the  vorticity  for  each  random 
mode. 


4.2.2  Spatial  Discretization 


Spatial  discretization  can  be  carried  out  by  any  method,  but  here  we  employ  the  spectral/hp  element  method  in 
order  to  have  better  control  of  the  numerical  error  [29].  In  addition,  the  all-spectral  discretization  in  space  and 
along  the  random  direction  leads  to  homogeneous  inner  products,  which  in  turn  results  in  more  efficient  ways  of 
inverting  the  algebraic  systems.  In  particular,  the  spatial  discretization  is  based  on  Jacobi  polynomials  on  triangles 
or  quadrilaterals  in  two-dimensions,  and  tetrahedra,  hexahedra  or  prisms  in  three-dimensions. 

4.3  Post-Processing 

The  coefficients  in  the  expansion  of  the  solution  process  (equation  (28))  are  obtained  after  solving  equations  (39) 
to  (43).  We  then  obtain  the  analytical  form  (in  random  space)  of  the  solution  process.  It  is  possible  to  perform  a 
number  of  analytical  operations  on  the  stochastic  solution  in  order  to  carry  out  other  analysis  such  as  the  sensitivity 
analysis.  Specifically,  the  mean  solution  is  contained  in  the  expansion  term  with  index  of  zero.  The  second-moment, 
i.e.,  the  covariance  function  is  given  by 


i?uu(xi,ti;x2,t2)  =  <  u(xi,ii)  -  u(xi,ti),u(x2,i2)  -  u(x2,t2)  > 

p 

=  'y]  [uj(xi, fi)uj(x2, f2)  <  4>2  >]  . 

i= 1 


(44) 
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Note  that  the  summation  starts  from  index  (i  =  1)  instead  of  0  to  exclude  the  mean,  and  that  the  orthogonality 
of  the  Askey-Chaos  basis  {<!>,;}  has  been  used  in  deriving  the  above  equation.  The  variance  of  the  solution,  i.e.  the 
‘mean-square’  value,  is  obtained  as 

Var{ u(x,  t))  =<  (u(x, t)  -  u(x,  t)J  >=  [uf (x, t)  <$?>],  (45) 

i=  1 

and  the  root-mean-square  ( rms )  is  simply  the  square  root  of  the  variance.  Similar  expressions  can  be  obtained  for 
the  pressure  field. 

5  Numerical  Results 


In  this  section  we  present  numerical  results  of  the  applications  of  the  generalized  polynomial  chaos.  We  first  consider 
the  stochastic  ordinary  differential  equation  and  demonstrate  exponential  convergence  with  the  optimal  Askey-Chaos. 
We  then  solve  the  incompressible  flow  in  the  pressure-driven  channel  where  there  is  uncertainty  associated  with  the 
wall  boundary  conditions.  Subsequently,  we  simulate  laminar  flow  past  a  circular  cylinder  with  uncertain  freestream. 

5.1  Stochastic  Ordinary  Differential  Equation 

5.1.1  Solution  Procedure 


To  demonstrate  the  convergence  type,  we  consider  the  ordinary  differential  equation 

dy(t) 


dt 


=  ~ky,  2/(0)  =  y , 


(46) 


where  the  decay  rate  coefficient  k  is  considered  to  be  a  random  variable  k{6)  with  certain  distribution  and  zero  mean 
value  k  =  0.  The  probability  function  is  f(k)  for  the  continuous  case  or  /(&*)  for  the  discrete  case.  The  deterministic 
solution  is  constant  over  time  y(t)  =  ye~kt  =  y,  while  the  mean  of  stochastic  solution  is 


y(t)=y  [  e  kt  f{k)dk  or  2/(t)=2/V/ 

Js 


,  —  kit 


m) 


(47) 


corresponding  to  the  continuous  and  discrete  distributions,  respectively.  The  integration  and  summation  are  taken 
within  the  support  of  the  corresponding  distribution,  and  in  general  the  mean  of  stochastic  solution  is  time-varying. 
By  applying  the  generalized  polynomial  chaos  expansion  of  equation  (19)  to  the  solution  y  and  random  input  k 


y{t)  =  ^2  k  =  ^2 


i= 0  i= 0 

and  substituting  the  expansions  into  the  governing  equation,  we  obtain 


Ef*--EE  >I>, 


2—0 


(48) 


(49) 


i— 0  j= 0 
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A  Galerkin  projection  onto  each  polynomial  basis  results  in  a  set  of  coupled  ordinary  differential  equations  for  each 
random  mode: 


dyijt) 

dt 


1 


p  p 


< 


<^2  >  y^,  eijihyj  (t)  i 


1  =  0,1,. 


.  p. 


(50) 


"  ^  0  j— 0 

where  =<  >.  A  standard  second-order  Runge-Kutta  scheme  is  used  to  integrate  the  equations.  We 

define  the  two  error  measures  for  the  mean  and  variance  of  the  solution 


i(*)  = 


V  (d)  2/exact  (t) 


r(*)  = 


tW 


'  exact 


is  the  variance  of  the  solution.  The 


2/exact  (f) 

where  y(t )  =  F%(f)]  is  the  mean  value  of  2/(t)  and  <r2(t)  =  |^(j/(f)  —  y(t)) 

initial  condition  is  fixed  to  be  y  =  1  and  the  integration  is  performed  up  to  t  =  1  (non-dimensional  time  units). 

5.1.2  Gaussian  Distribution  and  Hermite-Chaos 


(51) 


When  k  is  a  Gaussian  random  variable  with  probability  density  function  f(k)  =  /2,  the  optimal  Askey- 

Chaos  is  the  Hermite-Chaos  which  can  represent  the  input  k  ‘exactly’  with  first-order  expansion.  Figure  2  shows  the 
solution  by  the  Hermite-Chaos  expansion.  The  convergence  to  zero  of  errors  in  the  mean  and  variance  as  the  order 
of  Hermite-Chaos  increases  is  shown  on  semi-log  plot;  exponential  convergence  rate  is  achieved. 


Figure  2:  Solution  with  Gaussian  random  input  by  fourth-order  Hermite-Chaos;  Left:  Solution  of  each  random  mode, 
Right:  Error  convergence  of  the  mean  and  the  variance. 


5.1.3  Poisson  Distribution  and  Charlier-Chaos 

As  an  example  of  the  discrete  case,  we  assume  k  is  a  random  variable  with  Poisson  distribution 

\  ^ 

f(k;  A)  =  e-A  — ,  k  =  0, 1, 2, . . . ,  A  >  0.  (52) 

In  this  case  the  optimal  Askey-Chaos  is  the  Charlier-Chaos  (see  table  1).  Results  with  fourth-order  Charlier- 
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Figure  3:  Solution  with  Poisson  random  input  by  4th-order  Charlier-Chaos;  Left:  Solution  of  each  mode  (A  =  1), 
Right:  Error  convergence  of  the  mean  and  the  variance  with  different  A. 


Chaos  expansion  are  shown  in  figure  3  for  two  different  distributions  corresponding  to  different  values  of  A.  Again, 
exponential  convergence  rate  is  observed. 

5.1.4  Effects  of  Non-optimal  Basis 


In  this  section  we  present  examples  of  representing  a  stochastic  input  with  non-optimal  Askey-Chaos.  More  specif¬ 
ically,  we  present  results  of  using  Hermite-Chaos  expansion  for  an  exponential  distribution.  Although,  in  theory, 
Hermite-Chaos  converges  and  it  has  been  successfully  applied  to  some  non-Gaussian  processes  (e.g.,  lognormal  [24]), 
we  demonstrate  numerically  here  that  exponential  convergence  rate  is  not  realized.  In  figure  4  the  approximation  of 


p 


Figure  4:  Approximation  of  exponential  distribution  with  Hermite-Chaos;  Left:  PDF  of  different  orders  of  approx¬ 
imations  of  exponential  random  variable  by  Hermite-Chaos;  Right:  Error  convergence  of  the  mean  solution  with 
Laguerre-Chaos  and  Hermite  Chaos. 


an  exponential  random  variable  by  Hermite-Chaos  is  plotted  on  the  left.  It  can  be  seen  the  Hermite-Chaos  converges 
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and  the  fifth-order  approximation  is  very  close  to  the  exact  distribution,  with  some  noticeable  difference  at  x  ~  0 
where  the  PDF  reaches  its  peak  at  1.  Subsequently,  if  we  continue  to  use  Hermite-Chaos  to  solve  the  equation  (46) 
with  k  being  an  exponential  random  variable,  the  exponential  convergence  rate  will  not  be  maintained  as  opposed 
to  the  Laguerre-Chaos. 

Another  interesting  example  is  shown  in  figure  5,  when  a  beta  random  variable  is  approximated  by  the  Hermite- 
Chaos.  The  convergence  of  Hermite-Chaos  can  be  clearly  seen  from  the  approximated  PDF  compared  with  the 


Figure  5:  PDF  of  approximations  of  Beta  distributions  by  Hermite-Chaos;  Left:  a  =  /3  =  0,  the  uniform  distribution, 
Right:  a  =  2,  /?  =  0. 


exact  PDF.  The  special  case  of  a  =  j3  =  0  corresponds  to  the  uniform  distribution,  and  we  observe  oscillations 
near  the  corners  of  the  square.  This  is  analogous  to  the  familiar  Gibb’s  phenomenon  in  the  deterministic  spectral 
approximation.  In  this  case,  the  best  choice  is  the  Jacobi-Chaos  which  can  represent  the  beta  random  variable 
exactly  with  only  the  first-order  term.  We  expect  that  exponential  convergence  rate  will  not  be  maintained  if  the 
non-optimal  Hermite-Chaos  is  used  to  solve  equation  (46)  instead  of  the  Jacobi-Chaos. 

5.2  Pressure-Driven  Channel  Flow 

We  consider  a  pressure-driven  channel  flow  as  shown  in  figure  6,  where  the  boundary  conditions  are  considered  to  be 
uncertain.  The  domain  (see  figure  6)  has  dimensions  such  that  y  £  [—1,1]  and  x  £  [—5,5].  The  pressure  gradient, 
acting  like  a  driving  force,  is  equal  to  twice  the  kinematic  viscosity,  and  thus  for  a  no-slip  wall  condition  the  solution 
is  a  parabolic  profile  with  centerline  velocity  equals  unity. 

5.2.1  Pressure-Driven  Channel  Flow:  Uniform  Boundary  Conditions 

We  assume  that  the  boundary  conditions  at  the  two  walls  are  uncertain  with  zero  mean  value,  i.e. ,  u\  =  0  +  <ti£i  and 
U2  =  0  T  tj 2^2.  where  and  £2  are  two  idependent  random  variables,  and  and  02  are  their  corresponding  standard 
deviations.  Since  the  boundary  conditions  are  uniform  in  space,  with  periodic  boundary  conditions  specified  in  the 


17 


u=u2 


Figure  6:  Schematic  of  the  domain  for  pressure-driven  channel  flow  with  random  boundary  conditions. 

streamwise  direction,  the  nonlinear  terms  in  the  stochastic  Navier-Stokes  equations  (33)  vanish,  and  we  obtain  the 
exact  solution 

u(x,y)  =  (1  -  y2)  +  v{x,y)  =  0.  (53) 

The  solution  consists  of  a  parabolic  profile  for  the  mean  solution  and  two  linear  random  modes  (£i  and  <f2)  linearly 
distributed  across  the  channel  width.  Note  the  form  of  the  exact  solution  is  independent  of  the  distribution  type  of 
random  variables  and  £2. 


u0 

ul 


Figure  7:  Solution  of  the  pressure-driven  channel  with  uniform  Gaussian  random  boundary  conditions;  Left:  the 
solution  profile,  Right:  development  of  random  modes  of  v- velocity  with  nonzero  initial  conditions. 


On  the  left  of  figure  7  we  show  the  solution  profile  across  the  channel.  The  and  £2  are  two  independent 
Gaussian  random  variables  with  er \  =  0.02  and  cr2  =  0.01.  The  two-dimensional  (n  =  2)  Hermite-Chaos,  the  optimal 
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Askey-Chaos  in  this  case,  is  employed.  Although  the  solution  suggests  that  only  a  first-order  expansion  ( p  =  1)  is 
needed,  higher-order  terms  ( p  >  1)  are  included  in  the  computation  but  are  identically  zero  as  expected.  Another 
test  is  to  set  the  initial  condition  of  the  flow  to  an  arbitrary  random  state.  We  add  perturbation  terms  to  the  exact 
solution  (equation  (53))  for  each  random  mode  in  the  form  of  Uk(x,y,0)  =  apf(x,y)  and  Vk(x,y,  0)  =  apg(x,y)  for 
k  =  0, . . . ,  P.  Here  p  is  the  order  of  the  chaos  expansions  and  0  <  a  <  1  to  ensure  the  decaying  of  the  perturbation. 
On  the  right  of  figure  7  we  show  the  time  history  of  some  dominant  random  modes  of  u-velocity  at  the  center  of  the 
channel.  It  is  seen  that  due  to  the  nonlinear  interactions  between  the  random  modes  some  of  them  are  amplified  in 
the  early  stage,  but  eventually  all  modes  converge  to  the  exact  solution. 

Computations  with  other  types  of  random  inputs  have  been  conducted  with  their  corresponding  Askey-Chaos 
expansions.  More  specificly,  we  set  and  £2  to  be  beta  and  gamma  random  variables  and  employ  the  Jacobi-Chaos 
and  Laguerre-Chaos,  respectively.  Similar  results  were  obtained  with  the  results  shown  in  figure  7. 

5.2.2  Pressure-Driven  Channel  Flow:  Non-uniform  Boundary  Conditions 


Next  we  consider  the  case  of  non-uniform  random  boundary  conditions,  i.e.  the  wall  boundary  conditions  at  different 


locations  are  partially-correlated.  The  wall  boundary  conditions  are  assumed  to  be  random  processes  with  correlation 


function  in  the  form 


C(x  1,  £2)  =  u2e 


(54) 


where  b  is  the  correlation  length. 


This  correlation  function  has  been  employed  extensively  to  model  processes  in 


many  fields;  it  is  employed  here  because  it  allows  us  to  solve  the  eigenvalue  problem  (24)  of  the  Karhunen-Loeve 


expansion  of  equation  (23)  analytically.  If  this  is  not  the  case,  a  standard  numerical  eigenvalue  solver  can  be  used. 
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Figure  8:  Deviation  of  mean  solution  from  a  parabolic  profile  in  pressure-driven  channel  flow  with  partially-correlated 
random  boundary  conditions  at  the  lower  wall;  Upper:  u-velocity,  Lower:  u-velocity. 
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By  setting  a  relatively  large  correlation  length  b  =  100,  the  eigenvalues  of  the  Karhunen-Loeve  expansion  are 
Ai  =  9.675354,  A2  =  0.1946362,  A3  =  0.05014117, 

Due  to  the  fast  decay  of  the  eigenvalues,  we  use  the  first  two  terms  in  the  Karhunen-Loeve  expansion  given  by 
equation  (23).  This  results  in  a  two-dimensional  chaos  expansion  (n  =  2).  Resolution-independence  checks  were 
conducted  and  the  fourth-order  chaos  expansion  ( p  =  4)  were  found  to  be  sufficient  to  resolve  the  problem  in  the 
random  space.  Using  equation  (29)  this  results  in  a  fifteen-term  expansion  (P  =  14).  Only  the  lower  wall  boundary 
condition  is  assumed  to  be  uncertain  with  a  =  0.1,  while  the  upper  wall  is  stationary  and  deterministic.  A  parabolic 
velocity  profile  is  specified  at  the  inlet  and  zero  Neumann  condition  at  the  outlet.  A  mesh  with  10  x  2  elements  is 
employed  and  basis  Jacobi  polynomials  of  sixth-order  in  each  element  results  in  resolution  independent  solution  in 
space. 


Figure  9:  Contours  of  rms  of  it-velocity  (upper)  and  u-velocity  (lower). 


We  first  consider  the  lower  wall  boundary  condition  a  Gaussian  random  process  and  employ  the  Hermite-Chaos 
expansion.  Figure  8  shows  the  velocity  contour  plot  of  the  deviation  of  the  mean  solution  at  steady-state  from  a 
parabolic  profile.  The  mean  of  it-velocity  remains  close  to  the  parabolic  shape  and  the  mean  of  v- velocity,  although 
small  in  magnitude,  is  non-zero.  Figure  9  shows  steady-state  solutions  of  the  rms  (root-mean-square)  of  u  and 
u-velocity.  We  see  the  development  of  a  ‘stochastic  boundary  layer’  close  to  the  lower  wall.  All  the  higher-order 
expansion  terms  are  non-zero,  which  implies  that  although  the  random  input  is  a  Gaussian  process,  the  solution 
output  is  not  Gaussian.  Since  no  analytic  solution  is  available,  Monte  Carlo  (MC)  simulation  is  used  to  validate  the 
result.  Figure  10  shows  the  solution  of  mean  velocity  u  and  v  along  the  centerline  of  the  channel.  It  is  seen  that 
the  Monte  Carlo  solution  converges  non-monotonically  to  the  Hermite-Chaos  result  as  the  number  of  realizations 
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Figure  10:  Monte  Carlo  (MC)  and  Hermite-Chaos  (HC)  solution  of  the  mean  velocities  along  the  centerline  of  the 
channel;  Left:  w-velocity,  Right:  v-velocity. 

increases.  In  this  case,  it  is  only  after  40, 000  realizations  that  Monte  Carlo  solution  can  capture  the  solution 
accurately,  especially  the  nonlinear  interactions  close  to  the  inlet.  The  polynomial  chaos  solver,  with  15  terms  in  the 
expansions,  is  more  than  two  thousands  times  faster  than  the  Monte  Carlo  computation  without  using  any  special 
optimization  techniques.  In  figure  11  the  solution  of  the  mean  velocity  along  the  centerline  is  shown  corresponding  to 


Figure  11:  Hermite-Chaos  solution  of  the  mean  velocities  along  the  centerline  of  the  channel  with  different  cr;  Left: 
u- velocity,  Right:  v- velocity. 


different  values  of  cr.  It  can  be  seen  that  as  the  intensity  of  the  input  uncertainty  cr  increases  the  stochastic  solution 


responses  increase  nonlinearly. 


In  figure  12  we  plot  the  mean  solution  along  the  centerline  of  the  channel  with  different  types  of  stochastic  inputs. 


Specifically,  we  assume  the  random  processes  of  the  low  wall  boundary  condition  are  zero-mean  Gaussian,  uniform  and 
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exponential  processes  with  the  same  exponential  correlation  structure  (equation  (54))  and  fixed  parameter  a  =  0.4. 
The  corresponding  Askey-Chaos,  i.e. ,  the  Hermite-,  Legendre-  and  Laguerre-Chaos,  respectively,  are  employed.  The 
variance  of  the  velocity,  non-dimensionalized  by  the  input  variance  cr2,  is  shown  in  figure  13.  It  is  seen  that  the 
uniform  random  process  results  in  a  smoother  solution  with  smaller  variances  due  to  the  fact  that  the  uniform 
distribution  has  finite  support. 


X 


Figure  12:  Chaos  solution  of  mean  velocities  along  the  centerline  of  the  channel  with  different  types  of  input  processes; 
Left:  u- velocity,  Right:  v- velocity. 
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Figure  13:  Chaos  solution  of  variance  along  the  centerline  of  the  channel  with  different  types  of  input  processes;  Left: 
variance  of  w-velocity,  Right:  variance  of  w-velocity. 


Figure  14  shows  the  solution  of  mean  velocity  along  the  centerline  of  the  channel  corresponding  to  uniform 
stochastic  process  as  the  lower  wall  boundary  conditions,  with  the  same  correlation  structure  as  above  (cr  =  0.4). 
The  Legendre-Chaos  expansion  is  employed.  The  Monte  Carlo  solution  converges  to  the  chaos  solution;  with  120, 000 
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realizations  it  captures  the  nonlinear  interactions  near  the  inlet  accurately.  The  Legendre-Chaos  corresponds  to 
dimension  n  =  2  and  polynomial  order  p  =  4,  which  according  to  the  formula  of  equation  (29)  gives  15  terms  in  the 
expansion. 


Figure  14:  Monte  Carlo  (MC)  and  Legendre-Chaos  solution  of  the  mean  velocities  along  the  centerline  of  the  channel 
with  uniform  stochastic  inputs;  Left:  it-velocity,  Right:  ^-velocity. 


5.3  Flow  Past  a  Circular  Cylinder 

In  this  section  we  simulate  two-dimensional  incompressible  flow  past  a  circular  cylinder  with  random  fluctuations 
superimposed  to  the  freestream.  More  specifically,  the  inflow  takes  the  form  mn  =  u+g ,  where  g  is  a  random  variable 
or  process.  Here  we  focus  on  the  Gaussian  process  and  Hermite-Chaos  solution.  The  size  of  the  computational  domain 
is  [—15,  25]  x  [—9,  9]  and  the  cylinder  is  at  the  origin  (0,  0)  with  diameter  D  =  1.  The  definition  of  Reynolds  number 
is  based  on  the  mean  value  of  the  inflow  velocity  u  and  the  diameter  of  the  cylinder.  The  domain  consists  of  412 
triangular  elements  with  periodic  conditions  specified  in  the  crossflow  direction.  Sixth-order  Jacobi  polynomial  in 
each  element  is  observed  to  result  in  resolution-independent  solution  in  space  for  Reynolds  number  less  than  200. 

5.3.1  Flow  Close  to  First  Critical  Reynolds  Number 

It  is  well  known  that  for  two-dimensional  flow  past  a  circular  cylinder,  the  first  critical  Reynolds  number  is  around 
Re  ~  40,  where  the  flow  bifurcates  from  steady  state  to  periodic  vortex  shedding  [30].  Here  we  study  the  effects  of 
the  upstream  random  perturbations  close  to  this  Reynolds  number.  We  set  uln  =  u  +  cr£,  where  £  is  a  Gaussian 
random  variable  and  cr  is  its  standard  deviation.  The  one-dimensional  Hermite-Chaos  expansion  is  thus  employed. 
The  pressure  at  the  rear  stagnation  point  of  the  cylinder  is  extremely  sensitive  to  the  vortex  shedding  state  and  is 
monitored  in  our  computation. 
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Figure  15:  Time  history  of  mean  pressure  at  the  rear  stagnation  point  at  Re  =  40  (Gaussian  perturbation  with 
cr  =  0.1);  Left:  The  time  history,  Right:  Close-up  view. 

Figure  15  shows  the  time  history  of  the  mean  pressure  at  the  rear  stagnation  point  at  Re  =  40,  which  is  close 
to  the  critical  Reynolds  number.  Solution  with  fourth-order  and  sixth-order  Hermite-Chaos  are  shown,  together 
with  the  deterministic  pressure  history  as  reference.  A  negligible  difference  is  observed  between  fourth-order  and 
sixth-order  chaos  solutions  (less  than  0.1%).  Thus,  the  solution  can  be  considered  as  resolution-independent  in  the 
random  space.  In  the  close-up  view  we  see  that  the  10%  random  perturbation  (a  =  0.1)  triggers  an  instability  and 
the  flow  becomes  weakly  periodic,  as  opposed  to  the  deterministic  solution  which  remains  steady. 

Next,  we  lower  further  the  inflow  Reynolds  number  to  Re  =  35.  In  figure  16  we  show  the  time  history  of  the 
mean  pressure  signal  at  the  rear  stagnation  point.  Again,  resolution  independence  checks  show  a  negligible  difference 
(less  than  0.1%)  in  the  solutions  by  fourth-order  and  sixth-order  Hermite-Chaos.  It  is  shown  that  at  this  Reynolds 
number  a  10%  random  perturbation  (cr  =  0.1)  is  unable  to  trigger  an  instability  and  the  flow  remains  steady.  On 
the  other  hand,  with  a  larger  perturbation  (cr  =  0.2)  the  flow  becomes  weakly  unsteady  again. 

These  results  suggest  that  the  inflow  random  perturbations  have  noticeable  effects  on  the  stability  of  the  flow  near 
its  critical  Reynolds  number.  In  fact,  the  existence  of  upstream  perturbation  induces  the  instability  and  forces  the 
transition  to  occur  at  lower  Reynolds  number.  This  study  is  similar  to  that  of  [31]  where  the  convective  instability 
is  studied  by  introducing  random  perturbations  at  the  inflow  of  the  backward-facing  step  flow.  Instead  of  running 
many  realizations  of  the  deterministic  flow  solver,  here  we  can  resolve  the  propagation  of  inflow  uncertainty  by  chaos 
expansion  in  one  single  run  of  the  stochastic  solver. 
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Figure  16:  Time  history  of  mean  pressure  at  the  rear  stagnation  point  at  Re  =  35;  Left:  The  time  history,  Right: 
Close-up  view. 

5.3.2  Flow  at  Re  =  100 

We  consider  another  case  at  Re  =  100  with  freestream  random  velocity  partially  correlated.  The  inflow  is  U;vn  = 
u+g(y)  where  g(y)  is  a  Gaussian  process  with  the  exponential  covariance  kernel  of  equation  (54)  with  a  =  0.02.  Again, 
a  relatively  large  correlation  length  is  chosen  ( b  =  100)  so  that  the  first  two  eigenmodes  are  adequate  to  represent  the 
process  by  Karhunen-Loeve  expansion  (23).  Thus,  we  employ  a  two-dimensional  Hermite-Chaos  expansion  ( n  =  2) 
and  fourth-order  chaos  (p  =  4). 

Figure  17  shows  the  pressure  signal,  together  with  the  deterministic  signal  for  reference  (denoted  as  Pd  in  dotted 
line).  We  see  that  the  stochastic  mean  pressure  signal  has  a  smaller  amplitude  and  is  out-of-phase  with  respect  to 
the  deterministic  signal.  Although  initially,  the  stochastic  response  follows  the  deterministic  response,  eventually 
there  is  a  change  in  the  Strouhal  frequency  as  shown  in  figure  18.  Specifically,  the  Strouhal  frequency  of  the  mean 
stochastic  solution  is  slightly  lower  than  the  deterministic  one  and  has  a  broader  support. 

In  figure  19  we  present  velocity  profiles  along  the  centerline  for  the  deterministic  and  the  mean  stochastic  solution 
at  the  same  time  instant.  We  see  that  significant  quantitative  differences  emerge  even  with  a  relatively  small  2% 
uncertainty  in  the  freestream.  In  figure  20  we  plot  instantaneous  vorticity  contours  for  the  mean  of  the  vorticity  and 
compared  it  with  the  corresponding  plot  from  the  deterministic  simulation;  we  observe  a  diffusive  effect  induced  by 
the  randomness.  In  figure  21  we  plot  contours  of  the  corresponding  rms  of  vorticity.  It  shows  that  the  uncertainty 
influences  the  most  interesting  region  of  the  flow,  i.e. ,  the  shear  layers  and  the  near-wake  but  not  the  far-held. 
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Figure  17:  Pressure  signal  of  cylinder  flow  with  non-uniform  Gaussian  random  inflow.  Upper:  High  modes.  Lower 
Zero  mode  (mean). 


Figure  18:  Frequency  spectrum  for  the  deterministic  (high  peak)  and  stochastic  simulation  (low  peak). 


Figure  19:  Instantaneous  profiles  of  the  two  velocity  components  along  the  centerline  (in  the  wake)  for  the  determin¬ 
istic  and  the  mean  stochastic  solution. 


Figure  20:  Instantaneous  vorticity  field  :  Upper  -  Deterministic  solution  with  uniform  inflow;  Lower  -  Mean  solution 
with  non-uniform  Gaussian  random  inflow. 


Figure  21:  Instantaneous  contours  of  rins  of  vorticity  field  with  non-uniform  Gaussian  random  inflow. 

6  Summary  and  Discussion 

We  have  developed  a  stochastic  spectral  method  to  model  uncertainty  and  its  propagation  in  simulations  of  incom¬ 
pressible  flows.  Numerical  examples  were  presented  for  uncertain  boundary  conditions  but  the  method  can  also  be 
applied  to  model  uncertainty  in  the  boundary  domain,  e.g.  a  rough  surface,  in  the  transport  coefficients,  e.g.  the 
eddy  viscosity  in  large  eddy  simulations  and  other  transport  models,  or  in  interaction  forces  for  coupled  problems. 
It  provides  a  formal  procedure  for  constructing  a  composite  error  bar  for  CFD  applications,  as  proposed  in  [32] ,  that 
includes,  in  addition  to  the  discretization  errors,  contributions  due  to  imprecise  physical  inputs  to  the  simulation. 

More  specifically,  we  have  generalized  the  original  polynomial  chaos  idea  of  Wiener  and  proposed  a  broader 
framework,  i.e.  the  Askey-Chaos,  which  includes  Wiener’s  Hermite-Chaos  as  a  subset.  For  the  most  commonly 
known  probability  distributions,  there  exists  a  corresponding  “best”  orthogonal  functional,  in  the  sense  that  it 
leads  to  the  substantial  dimensional  reducibility.  We  also  applied  the  Askey-Chaos  expansion  to  the  Navier-Stokes 
equations  to  model  uncertainty  in  incompressible  flow  simulations  for  steady  and  unsteady  problems.  Convergence 
was  verified  with  comparisons  against  exact  solutions  and  solutions  from  Monte  Carlo  simulations  for  steady-state 
problems. 

As  regards  efficiency,  a  single  Askey-Chaos  based  simulation,  albeit  computationally  more  expensive  than  the 
deterministic  Navier-Stokes  solver,  is  able  to  generate  the  solution  statistics  in  a  single  run.  Specifically,  an  Askey- 
Chaos  Navier-Stokes  simulation  with  a  total  number  of  terms  (P  + 1)  is  approximately  (P  +  1)  times  more  expensive 
than  the  corresponding  deterministic  one.  (The  overhead  associated  with  the  coupling  terms  is  negligible).  In 
contrast,  for  the  Monte  Carlo  simulation,  tens  of  thousands  of  realizations  are  required  for  converged  statistics,  which 
is  prohibitively  expensive  for  most  CFD  problems  in  practice.  For  the  problems  we  studied  here  we  can  only  make 
direct  comparisons  for  the  steady-state  cases,  e.g.  the  pressure-driven  flow  with  random  boundary  conditions.  In  the 
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case  of  the  Hermite-Chaos,  the  speed-up  factor  is  approximately  2,500  (see  figure  10;  (40,000/15))  for  comparable 
accuracy  of  the  first  two  moments.  Similarly,  for  the  Legenclre-Chaos  the  speed-up  factor  is  about  8,000  (see  figure 
14;  (120,000/15)).  Clearly,  better  Monte  Carlo  algorithms  (e.g.,  accelerated  versions  with  variance  reduction)  could 
reduce  this  factor  but  we  still  expect  a  three  orders  of  magnitude  speed-up  for  most  problems  involving  a  relatively 
low  number  of  random  dimensions. 

There  are,  however,  several  open  problems  that  need  to  be  resolved  for  the  Askey-Chaos  to  be  a  robust  and 
effective  tool  for  modeling  uncertainty.  In  particular,  further  research  is  required  on: 

•  Convergence  rate.  The  relatively  poor  resolution  properties  of  Hermite  and  Laguerre  expansions,  compared  to 
other  spectral  polynomials,  are  well  documented  in  the  literature  [33,  34].  However,  re-scaling  procedures,  as 
done  in  [35],  can  be  applied  or  a  change  of  the  trial  basis  from  the  Askey  scheme,  as  demonstrated  in  section 
5.1.4,  can  be  employed  to  accelerate  convergence. 

•  Dimensionality  of  the  stochastic  input.  This,  in  turn,  determines  the  dimensionality  of  the  random  space  and 
correspondingly  the  computational  complexity  of  the  problem.  For  a  physical  input  random  process  with  a 
very  short  correlation  length,  a  high  dimensional  chaos  expansion  is  required.  As  shown  in  equation  (29),  the 
number  of  expansion  terms  ( P  +  1)  increases  fast,  although  algebraically,  both  with  the  dimension  n  as  well 
as  the  polynomial  order  p.  In  contrast,  the  convergence  rate  of  the  Monte  Carlo  method  is  independent  of  the 
number  of  random  dimensions. 

•  Interaction  with  spatial/temporal  discretization.  In  this  paper,  we  have  employed  a  multi-step  time  integrator 
and  spectral//ip  element  methods  for  discretizing  the  deterministic  operators.  Some  of  the  conclusions  of  the 
work  presented  here  may  not  be  readily  extended  to  other  discretizations. 

The  aforementioned  issues  can  be  addressed  by  systematic  studies,  investigating,  for  example,  different  type 
projections,  filtering,  proper  trial  basis,  rescaling,  etc.  We  are  currently  working  in  addressing  some  of  these  issues 
and  we  will  report  results  in  future  publications. 
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A  Some  Orthogonal  Polynomials  in  Askey  scheme 


In  this  section  we  briefly  review  the  definitions  and  properties  of  some  important  orthogonal  polynomials  from  Askey 
scheme,  which  are  discussed  in  this  paper  for  the  Wiener-Askey  polynomial  chaos. 

A.l  Continuous  Polynomials 

A. 1.1  Hermite  Polynomial  Hn(x)  and  Gaussian  Distribution 


Definition: 


Orthogonality: 


Recurrence  relation: 


Hn(x) 


(2x)n  2F0 


e  x  Hm(x)Hn(x)dx  =  2 nn\8„ 


Hn+ i(x)  -  2 xHn(x)  +  2nHn_i(x)  =  0. 


(55) 


(56) 


(57) 


Rodriguez  formula: 

e-**Hn{x)  =  {-ir^(e-*2).  (58) 

The  weighting  function  is  w(x)  =  e~x  from  the  orthogonality  condition  (56).  After  rescaling  x  by  y/2,  the 
weighting  function  is  the  same  as  the  probability  density  function  of  a  standard  Gaussian  random  variable  with  zero 
mean  and  unit  variance. 


A. 1.2 


Laguerre  Polynomial  L 


(a) 


{x) 


and  Gamma  Distribution 


Definition: 

L(n\x)  =  ^  iTi(— n;o  +  1;  x). 

n\ 

Orthogonality: 

[  e~xxaL^\x)L^\x)dx  =  T(n  +  °  +  ^  Smn,  a  >  -1. 
Jo  n\ 

Recurrence  relation: 

(n  +  1  )L{n+i(x)  -  (2n  +  a  +  1  -  x)L^\x)  +  (n  +  a)!^^)  =  0. 

Rodriguez  formula: 

1  dn 

e~xxaL^\x)  =——  (e~xxn+a)  . 
n  w  n\  dxn  v  ’ 


(59) 


(60) 


(61) 


(62) 
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Recall  that  the  Gamma  distribution  has  the  probability  density  function 


fix)  = 


o-x/P 


a  >  — 1,  /3  >  0. 


(63) 


/?“+1r(a  +  l)’ 

Despite  of  the  scale  parameter  /?  and  a  constant  factor  T(a+  1),  it  is  the  same  as  the  weighting  function  of  Laguerre 
polynomial. 


A. 1.3  Jacobi  Polynomial  and  Beta  Distribution 


Definition: 


O  r  t  hogonality : 


where 


2  n  ( 


— n,  n  a  (3  1;  a  1: 


1  —  x 


J  (l-a:)“(l  +  xf  P^’P\x)p(a’0Xx)dx  =  h2n8mn,  a>- 1,  /?  > -1, 

2«+/3+i  r(n  +  a  +  l)r(n  +  /3  +  1) 


= 


2?r  +  a  +  /3  +  1  T(n  +  a  +  /3  +  l)n! 


Recurrence  relation: 


xP^\x)  = 


2 {n  l)(n  +  (x  -\-  [3  1) 

(2n  +  a  +  (3  +  l)(2n  +  a  +  f3  +  2) J  n+1 
/32  —  a2 


(2t&  -b  ck  -b  /3)(277-  -b  u  -b  f3  -b  2) 

2(n  +  a)(n  +  /3)  p(a,/3) 

(2t&  -b  q:  ~b  /3)(2?7-  +  a;  +  /?  +  1) 


p:p\x) 

d-f’w- 


Rodriguez  formula: 


rln 

(l  -  *)“(i  +  xfP^Xx)  =  ^  —  [(l  -  z)n+“(i  +  *r+/3] . 


2nn!  da;” 

The  Beta  distribution  has  the  probability  density  function 

( x  —  a)P{b  —  x)“ 


/(z)  = 


(b  —  a)a+P+1B(a  +  1,  /?  +  1) : 
where  B(p,q)  is  the  Beta  function  defined  as 

s(p,  S)  =  r(p)r(«) 


a  <  x  <  b, 


(64) 


(65) 


(66) 


(67) 


(68) 


(69) 


r  ip+q)  ' 

It  is  clear  that  despite  of  a  constant  factor  the  weighting  function  of  Jacobi  polynomial  w(x)  =  (1  —  x)“(l  +  x)@  from 
(65)  is  the  same  as  the  probability  density  function  of  Beta  distribution  defined  in  domain  [—1, 1],  When  a  =  (3  =  0, 
the  Jacobi  polynomials  become  the  Legendre  polynomials  and  the  weighting  function  is  a  constant  which  corresponds 
to  the  important  uniform  distribution. 
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A. 2  Discrete  Polynomials 

A. 2.1  Charlier  Polynomial  Cn(x;a)  and  Poisson  Distribution 


Definition: 


Orthogonality: 


Recurrence  relation: 


Cn(x;a)  =  2F0  -n,-x\  . 

a  ' 


^  .  Cm{x ,  Cl^Cn(x,  tt)  —  CL  C  TL ^$rnm  Cl  0. 

Xl 

x—0 


—xCn(x\  a)  =  aCn+i(x;  a)  —  (n  +  a)Cn(x ;  a)  +  nCn_i(x]  a). 


Rodriguez  formula: 

where  V  is  the  backward  difference  operator  defined  as 

A/O)  =  f(x  +  1)  -  f(x)  and  V/O)  =  /O)  ~  /O  -  !)■ 


The  probability  function  of  Poisson  distribution  is 

/0;a)  =  e_o^T,  fc  =  0,l,2, .... 
xl 

Despite  of  a  constant  factor  e~a,  it  is  the  same  as  the  weighting  function  of  Charlier  polynomials. 

A. 2. 2  Krawtchouk  Polynomial  Kn(x;p,  N)  and  Binomial  Distribution 


(70) 

(71) 

(72) 

(73) 

(74) 

(75) 


Definition: 

Kn(x;p,N)  =  2Fi  (-n,  -x;  -N;  ^  ,  n  =  0,l,...,N. 

Orthogonality: 

X]  (^)pX(1  ~P)N~XKm(x;P,  N)Kn(x;p,  N)  =  (~T")  6rnn,  0  <P<  1- 

Recurrence  relation: 


(76) 


(77) 


- xK(x;p,N )  =  p(N  -  n)Kn+i(x;p,  N)  -  \p(N  -  n)  +  n(l  -p)\Kn(x;p,  N) 

+  n{l-p)Kn-i{x\p1N). 


(78) 
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Rodriguez  formula: 


N\  (  p 


x  )  V1  ~P 


Kn{x\  p,  N)  =  V‘ 


N  —  n\  /  p 
x  )  V1  -p 


Clearly,  the  weighting  function  from  (77)  is  the  probability  function  of  the  binomial  distribution. 

A. 2. 3  Meixner  Polynomial  Mn(x\f3,c)  and  Negative  Binomial  Distribution 


Definition: 


Orthogonality: 


Mn{x\p,c)  =  2-Fi  l-n,  /?;  1 - . 

V  c 


YJ  — -cxMm(x;  p,  c)Mn(x ;  p,  c)  = 


(/?)„(! -c)0 


&7nm  P  P  0;  0<C^1. 


Recurrence  relation: 


{c-l)xMn(x\P,c)  =  c(n  + P)Mn+i(x-,P,c)  ~[n+ (n  + P)c]Mn(x-,P,c) 
+  P,  c). 


Rodriguez  formula: 


The  weighting  function  is 


f{x)={-^f{l-cfcx,  0<p<l,  /?  >  0,  *  =  0,1,2,....  (84) 

x\ 

It  can  verified  that  it  is  the  probability  function  of  negative  binomial  distribution.  In  the  case  of  P  being  integer,  it 
is  often  called  the  Pascal  distribution. 

A. 2. 4  Hahn  Polynomial  Qn(x\a,  P,N)  and  Hypergeometric  Distribution 


Definition: 


Qn(x-,a,P,N)=  3F2(—n,n  +  a  +  P+l,—x-,a  +  l,—N-,l),  n  =  0,l,...,N. 


Orthogonality:  For  a  >  —  1  and  P  >  —  1  or  for  a  <  — N  and  P  <  — N , 


a  +  x\fP  +  N  —  x 


x  )  V  N  —  x 


Qm(x,  O,  /?,  -/V)t^n(x,  Q:,  /?,  A) 
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where 


2  _  (— l)"(?r  +  a  +  (3  +  1)n+i{(3  +  1  )n'n\ 
n  ~  (2 n  +  a  +  (3  +  1  )(a  +  l)„(-iV)„IV!  ' 

Recurrence  relation: 


xQn{x3)  —  ( d ri  -\-  Cji)Qn(x)  3~  CnQn—  l(x), 


(87) 


where 


Qn(x)  :=  Qn(x;a,f3,N) 


and 


Rodriguez  formula: 


_  (n+a+/3+l)(n+a+l)(N—n) 

n  ~  (2n+o:+/3+l)(2n+a+/3+2) 


_  rt(ra+a+,3+iV+l)(re+/3) 
—  (2n+a+/3)(2n+a+/3+l)  ' 


w(x;  a,  f3,  N)Qn(x\  a,  /?,  IV) 


(~ l)”(/3  +  l)r 

(-JV)n 


Vra[w(:r;  a  +  n,  (3  +  n,  N  —  n)], 


(88) 


where 


w(x\  a,  (3,  N)  = 

If  we  set  a  =  —a  —  1  and  (3  =  —  (3  —  1,  we  obtain 


a  +  af\  f  f3  +  N  —  x 
N-x 


){x)  = 


i  a(/-j 

r"V_i)  re5) 


Apart  from  the  constant  factor  l/(  this  is  the  definition  of  hyper  geometric  distribution. 
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